function [J1 J2 J1_1 J2_1] = func_J_k(nxy, nx, CONSTS)

    nz = func_nz_from_nxy(nxy, CONSTS);

    [nxy1 nxy2] = func_nxyk_from_nz(nz, CONSTS); 
    [ny1 ny2] = func_nyk_from_nxyk(nz, nx, CONSTS);
    [alfa_1 alfa_2 beta_1 beta_2] = func_alfa_beta_from_nz(nz, CONSTS);
        
    J1 = (1i*ny1+alfa_1.*nx)./(nxy1.^2);
    J2 = (1i*ny2+alfa_2.*nx)./(nxy2.^2);
    J1_1 = (1i*ny1-beta_1.*nx)./(nxy1.^2);
    J2_1 = (1i*ny2-beta_2.*nx)./(nxy2.^2);
    
end
